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The turbulent flows through a plane channel and a channel with a constriction (2-D 
hill) are numerically simulated using DNS and RANS calculations. The Navier-Stokes 
equations in the DNS are solved using a higher order kinetic energy preserving central 
schemes and a fifth order accurate upwind biased WENO scheme for the space 
discretization. RANS calculations are performed using the NASA code CFL3D with the k- 
omega SST two-equation model and a full Reynolds stress model. Using DNS, the 
magnitudes of different terms that appear in the enstrophy equation are evaluated. The 
results show that the dissipation and the diffusion terms reach large values at the wall. All 
the vortex stretching terms have similar magnitudes within the buffer region. Beyond that 
the triple correlation among the vorticity and strain rate fluctuations becomes the 
important kinematic term in the enstrophy equation. This term is balanced by the viscous 
dissipation. In the separated flow, the triple correlation term and the viscous dissipation 
term peak locally and balance each other near the separated shear layer region. These 
findings concur with the analysis of Tennekes and Lumley, confirming that the energy 
transfer terms associated with the small-scale dissipation and the fluctuations of the vortex 
stretching essentially cancel each other, leaving an equation for the dissipation that is 
governed by the large-scale motion. 


I. Introduction 

T he major objectives of any CFD code are to predict global aerodynamic quantities such as lift and drag as well 
as to predict local flow features such as skin friction, heat transfer, pressure oscillations, etc. accurately and 
efficiently for a wide range of problems. The flow physics are fundamentally governed by the unsteady three- 
dimensional Navier-Stokes (N-S) equations that state the three conservation laws for mass, momentum, and energy. 
Turbulent flow is characterized by the existence of a vast range of length and time scales ranging from the smallest, 
the Kolmogorov scale, to the largest, determined by the geometry. 1 ' 2 The required computer resources and time 
constraints hinder the solution of these equations numerically for turbulent flows at high Reynolds numbers. 
Analysis 2 has shown that the number of grid points required increases as Re 8 3 and the number of time steps 
required increases as Re 4 ' 3 . These challenging requirements restrict the solution of the full N-S equations to simple 
geometries such as channels and flat plates at low Reynolds numbers. The first successful direct numerical 
simulation'' (DNS) was performed for a channel flow at a Reynolds number of 3000. With increasing 
computational capability, the application of DNS to more complex problems at high Reynolds numbers is currently 
pursued by many researchers. The recent status of DNS in turbulent flow is reviewed by Moin and Mahesh. 4 
However, DNS usage is limited to understanding the turbulent physics in canonical problems and to extending the 
findings for applications to turbulence model development. Hence, the current state of the art is to solve some 
approximate versions of the N-S equations that require modeling. 

The simplest and the most popular approximate method is the set of Reynolds-Averaged Navier-Stokes 
(RANS) equations, where the equations are derived for time-averaged quantities. The unclosed Reynolds stresses 
are modeled to close the equations. Wilcox 5 gives a good account of different models and their applications. In 
general, existing models provide good results compared to experiments in attached flows where turbulence is in 
equilibrium. However, prediction of separated flows with existing RANS turbulence models is still not very 
satisfactory. It is difficult to accurately capture both the separation and reattachment points (i.e., the size of the 
separation bubble). In particular, it is believed that the flow field dynamics near the separation point are not 
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predicted correctly. The RANS computations performed for the flow over a two-dimensional hump, NASA 
Langley CFDVAL2004 Case 3 6 ' 7 , consistently over-predicted the reattachment point compared to the experiment. 
The computed Reynolds stresses were 2 to 3 times smaller than measured values. Computations performed 8 for two 
other massively-separated flow cases (the 2-D periodic hill 9 and the 3-D Ahmed body 10 ) using different turbulence 
models also predicted too-long separation bubbles. These computations further exemplify the deficiencies in 
predicting massively-separated flows using existing RANS models. 

The second approximate method uses the Large-Eddy-Simulation (LES) equations, which are obtained for the 
large-scale motions by filtering small scales. The unknown stresses introduced by the small scales (subgrid-scale 
stresses, or SGS) are modeled to close the equations. Sagaut 11 and Gamier et al. 12 give account of different SGS 
models and their applications in incompressible and compressible flows, respectively. This approach is 
computationally more expensive than the RANS approach, but yields more accurate results without many ad-hoc 
modifications to the SGS models. This methodology is being applied to a wide range of problems. Recently, 
several DNS and LES simulations have been performed for separated flows over humps in channel flows. 13 " 18 
Comparisons of the RANS solutions obtained using two-equation models and Reynolds-Stress models 19 with the 
LES and DNS solutions led to the same conclusion that the RANS solutions over-predict the separation bubble 
lengths. Another observation from the LES and DNS was that the ratio of production over dissipation takes very 
high values in the separated shear region compared to homogeneous flows. Whether this is one of the causes for the 
deficiency in RANS is not clear. The turbulent dynamics between the separation point and the reattachment point is 
very complex. 20 The adverse pressure gradient near the reattachment point and spreading of the shear layer due to 
turbulence are interconnected and the balance between these two factors determines the separation and the 
reattachment points. The focus of this paper is a numerical study of how turbulence production and dissipation 
behave in separated flows. 

The production is determined by the Reynolds stresses and the mean velocity gradients. It is generated by the 
large-scale dynamics of turbulence. The dissipation takes place at the smallest scales due to molecular viscosity. 
According to Kolmogorov’s hypothesis, the turbulent kinetic energy in the large scales is transferred to the small 
scales across the inertial range and is dissipated due to molecular viscosity in the smallest scales. This energy 
process, i.e., how the large scales and the small scales communicate, is an unsolved problem in turbulence. 
Tennekes and Lumley 21 described this energy transfer as a vortex stretching process due to nonlinearity. It is 
known that mean-square vorticity fluctuations (enstrophy) equal the energy dissipation in homogeneous flows. 
Hence an equation for the enstrophy involves the small scales as well as large-scale dynamics. It is postulated that 
the terms associated with the small-scale dissipation and the fluctuations of the vortex stretching are higher order in 
powers of the Reynolds number and will cancel each other, leaving an equation for the dissipation that is governed 
by the large-scale motion. One of our objectives is to verify this hypothesis numerically. The enstrophy equation 
consists of averages of second derivatives squared and second derivatives of higher order averages. They are 
difficult to compute accurately in high Reynolds number flows. 

As a first step, we performed a simulation for a periodic plane channel flow at a low Reynolds number and 
evaluated the terms in the enstrophy equation. The second step was to perform a simulation for a flow through a 
channel with a constriction (separated flow over a 2-D hill). This problem is a standard test case for the validation 
of codes and turbulence models. This case has been investigated numerically using LES 9 ’ 15 ’ 16,18 and DNS. 17 We also 
performed RANS calculations using two-equation and Reynolds stress models for the separated flows over the hill 
and compared the mean and the turbulence quantities with the DNS results. 

II. Governing Equations 

The partial differential equations solved are the three-dimensional unsteady compressible Navier-Stokes 
equations in conservation form 
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Here ( x,y,z ) are the Cartesian coordinates, (u,v,w) are the velocity components, p is the density, and p is the 
pressure. E is the total energy given by 

u 1 + V 2 + w 2 

E = e + , 

2 

e = cT, p =pRT. (3) 


Here e is the internal energy and T is the temperature. The shear stress and the heat flux are given by 
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The viscosity (p) is computed using the power law p = const. T°' 7 , and the coefficient of conductivity (k) is given 
in terms of the Prandtl number Pr. We performed the computations for a constant Prandtl number of 0.7. The flows 
are assumed to be periodic in the streamwise direction, hence the mean pressure gradient is zero. The flow is 
maintained by applying a body force fj(t) in the streamwise direction. This is a function of time only and is 
determined at every step of the time marching by requiring that the average mass-flux remains constant. 
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A. Solution Algorithm 

The governing equations were solved using two spatial discretization schemes. One is based on a 5 th -order 
accurate weighted essentially non-oscillatory (WENO) scheme. We have been using this scheme in high-speed 
boundary layer transition simulations. The WENO method and formulas are explained in Shu. 22 The application of 
the ENO method to the Navier-Stokes equations is presented in Atkins. 23 The solution method implemented in the 
present computations is described in Balakumar. 24 This scheme is based on an upwind biased algorithm and has 
more dissipation than central difference schemes in spite of having better stability characteristics than the central 
schemes. 

Recently, we added higher order kinetic energy preserving central schemes to the code. The methodology for 
the discretization of the Euler fluxes is the same as described in Ducros et al. 25 and Pirozzoli 26 The code can run to 
any order on uniform grids. The viscous terms have not been written in Laplacian form, and are currently solved 
using a fourth order scheme with double differentiation. When we tried the scheme in curvilinear grids as in the 2- 
D hill configuration, high frequency oscillations started to appear in the solution after a certain number of steps. 
Hence, we stopped using the central scheme for this case and performed the simulation with the higher order 
upwind scheme. The plane channel flow simulations were performed using an 8 th order central scheme in the 
streamwise and spanwise directions and a 4 th order scheme in the wall normal direction for the Euler fluxes and 4 th 
order central schemes for viscous fluxes. We employed a 3 rd -order total-variation-diminishing (TVD) Runge-Kutta 
scheme for time integration with both methods. 

We used a body-fitted curvilinear grid system in all of the simulations. The equations are transformed from the 
physical coordinate system (x, y, z) to the computational curvilinear coordinate system (£, q, Q. The grids were 
uniform in the streamwise and spanwise directions and were stretched in the normal directions close to the walls. 
We will provide more details in the following section. 
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III. Models and Flow Conditions 

Computations were performed for two flow geometries. First was a flow through a plane channel and second 
was a flow through a channel with a constriction (2D hill). The geometries and the coordinate systems are shown in 
Figs. 1(a) and 1(b). The flow was periodic in the streamwise and spanwise directions for the both cases. 


(a) Plane channel flow 

In the plane channel flow, the lengths are non-dimensionalized by the half width of the channel, h. The 
computational dimension in the streamwise direction is L x = 2n, in the spanwise direction is L z =4n/3, and in the 
normal direction is L y = 2. These are the same dimensions used in the DNS simulations 27 and in the LES 
simulations. 28 Simulations are performed for a Reynolds number based on the bulk velocity, density and the 
channel half width of Ret = pt Ut h/p w = 3000. The wall temperature is fixed at a constant temperature of 353K. 
The Mach number based on the bulk velocity is M 0 = 0.5. The non-dimensional bulk density and the bulk velocity 
are defined as 
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We used (1201, 251, 401) points in the streamwise, wall normal and spanwise directions, respectively. The grid 
spacings in viscous units are Ax + = 2 and Az + = 2 in the streamwise and spanwise directions, respectively. The 
minimum grid spacing at the wall is Ay + = 0.16 and the spacing at the center of the channel is Ay + = 3.4. 

(b) Plane channel flow with a constriction (2-D hill) 

In the 2D hill flow, the lengths are non-dimensionalized by the hill height, h. The computational dimension in 
the streamwise direction is L x = 9, in the spanwise direction is L z = 4.5, and in the normal direction above the crest 
of the hill is L y = 3.035. These are the same dimensions used in the DNS simulations 17 and in the LES simulations. 18 
Simulations are performed for a Reynolds number based on the bulk velocity at the entrance of the channel x = 0, 
density and the channel half width of Ret = pt Ut h/ u w = 2800. The wall temperature is fixed at a constant 
temperature of 353K. The Mach number based on the bulk velocity is M 0 = 0.2. The non-dimensional bulk density 
and the bulk velocity are defined as 



We used (801, 351, 513) points in the streamwise, wall normal and spanwise directions, respectively. The grid 
spacings in viscous units change along the streamwise direction due to the change in friction velocity. The maxium 
grid spacings in the x- and z-directions are Ax + = 4 and Az + = 3, respectively. The minimum grid spacing at the wall 
varies between Ay + = 0. 1 to 0.3 and the maximum spacing in the middle of the channel is Ay + = 4. 

IV. Results 


A. Plane channel 

(a) Mean quantities 

The statistical quantities are obtained by averaging the solution in the spanwise direction and in time. The 
averaging in time was performed for about 20 flow through periods. One flow through period is based on the 
channel length divided by the bulk velocity. Table I compares the global quantities with the DNS data of Kim et 
al. 3 and with that obtained using LES with the Smagorinsky modeF 8 for the same free stream Mach number of 0.5. 
The agreement is very good compared to other computations. The wall friction velocity is about 3% smaller than 
that obtained in the DNS of incompressible channel flow. 
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Table 1. Mean flow variables 



Present DNS 
M 0 = 0.5 

DNS 

Incompressible 
Ref. 3 

LES 

Smagorinsky 
Ref. 28 

u c 
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1.1630 

W, 

6. 19*1 O' 2 
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5.8521*10 2 

c / 

7.66*10" 2 

8.18*1 0" 3 

7. 1 454* 1 0" 3 

TJT w 

1.047 


1.047 

Pc ! Pw 
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Figures 2(a) and (b) compare the computed mean velocity profdes with the DNS data of Kim et al. 3 Figure 2(a) 
depicts the results non-dimensionalized by the bulk velocity and the channel half-width, while Fig. 2(b) shows the 
results non-dimensionalized by the viscous units. We also include the linear and the logarithmic profdes for 
comparison. The current computed results are about 2% higher than the logarithmic curve. We are pursuing another 
simulation with a finer grid distribution near the wall to check the grid sensitivity of the results. The half-width of 
the channel in y + units is 185.7. This length is small due to the low Reynolds number simulated in this case. The 
logarithmic layer extends from y + ~30 to 180. Figures 3(a) and (b) compare the turbulent stress intensities non- 
dimensionalized by the wall friction velocity with the DNS data of Kim et al. Fig. 3(a) shows the results for the 
intensities of the normal stresses, and Fig. 3(b) depicts the results for turbulent shear stress. The agreement is very 
good. Maximum intensities are slightly higher than those for the incompressible case. The ratios of the three 
normal stresses are u rms : v rms : w rms = 1: 1/3 :l./2.5. The maximum rms values occur at v + =15, 52 and 37 for the u, v 
and w components respectively. The maximum shear stress is 0.76 and this occurs aty + = 32. 

For a compressible flow, the transport equations for the Reynolds stresses and the turbulent kinetic energy take 
the following forms. These equations appear in several articles; we followed Canuto’s 29 derivation. 
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In these equations, an overbar denotes a Reynolds averaging and a tilde denotes a Favre averaging. Hence any 
variable q, is written as 


q = q+ q 
= q + q" 

The term on the left hand side of Eq. (8) represents the advection due to the mean flow. The terms on the right hand 
side are: (1) Sg. Production, (2) IIj, pressure-strain, (3) Dg, diffusion, (4) Sj, dissipation and (5) Bg, A, Fg, 
compressibility terms. The diffusion of Reynolds stress term, Dg, consists of three parts: (a) diffusion due to 
turbulence, (b) diffusion due to pressure, and (c) diffusion due to viscosity. There are three compressibility terms: 
two of them are associated with the Favre averaging term, and the other one is associated with the dilatation. The 
trace of this equation yields the equation for the turbulent kinetic energy: 
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The left hand side of Eq. (9) is the rate of change of turbulent kinetic energy and the terms on the right hand side 
are: (1) production, (2) diffusion, (3) dissipation, and (4) compressibility effects (B, A and F). The pressure strain 
term is a redistribution term for the Reynolds stresses and that term disappears in the kinetic energy equation. 
Figure 4 shows the magnitudes of the different terms in the kinetic energy equation for the channel flow. The 
results agree with the computational results of Moser et al. 27 Turbulent energy is transferred to the wall by diffusion 
and is dissipated by the molecular viscosity. The diffusion becomes negligible beyond about v + = 6.3 (y/h ~ 0.10) 
and the production and dissipation balance in this region. The production and diffusion peak at about v + =12 (y/h ~ 
0.10) and the magnitudes of the production, diffusion and dissipation at this location are 0.233, 0.104 and 0.128, 
respectively. The ratio of production to dissipation near the peak is about P/e =1.8. The “balance” represents the 
sum of all the individual terms, and is approximately zero as expected. The enstrophy (discussed in the next 
section) is also plotted in this figure. Figures 5(a-d) depict the balance of different terms in the Reynolds stress 
equations for uu, vv, ww and uv, respectively. The results agree with the computations of Mansour et al. 30 The 
production is positive only for the uu component and zero for the vv and ww components. The pressure-strain is 
opposite to this, being negative for the uu component and positive for the other two normal components. This is the 
well-known conclusion that the turbulence is produced by the mean strain in the streamwise direction and is 
transferred to the other two components by the turbulent pressure and strain rate fluctuations. 

(b) Enstrophy computations 

Vorticity fluctuations play a key role in transferring energy among different scales of the turbulent 
fluctuations. 21 Three processes that balance the turbulent dynamics are the vortex stretching, viscous diffusion and 
viscous dissipation. Vortex stretching is a kinematic process and can occur at the largest and the smallest scales of 
motion due to velocity strain. At high Reynolds numbers, viscous diffusion and viscous dissipation occur at small 
scales. Hence, vortex stretching is the coupling between the large scale and small-scale motions. Tennekes and 
Lumley 21 described the importance of the vortex stretching process for the transfer of energy from large scale to 
small scales by analyzing the order of magnitude of the terms that appear in the equation for the mean-square 
vorticity fluctuations, or the enstrophy. They showed that the stretching of vorticity by the turbulent fluctuations 
has the highest magnitude compared to other terms in the equation. This term is balanced by the viscous destruction 
term in the enstrophy equation. The enstrophy is equivalent to the turbulent kinetic energy dissipation in 
homogeneous flows, and these two terms are approximately equivalent in inhomogeneous flows. Hence our 
objectives are first to evaluate the terms in the equation for the enstrophy using the DNS data, and second to 
compare the model equation for the turbulent kinetic energy dissipation with the exact dissipation or the enstrophy 
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equation. Mansour et al. 30 performed a similar analysis for the plane channel flow using the DNS data of Kim et 
al. 3 The numerical calculations become difficult due to the appearance of higher order derivatives in the averaged 
terms for the enstrophy equation. We will first present the results for the channel flow from the present simulation 
data, and in the next section will present the results for the separated flow over a hill. The vorticity fluctuations 
amplitude is defined as the enstrophy: 


Enstrophy = 

The enstrophy is related to the dissipation term as (Morinishi et al. 31 ) 
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The enstrophy term is plotted in Fig. 4 for comparison with the total kinetic energy dissipation term e. The figure 
shows that the enstrophy and the turbulent kinetic energy dissipation are almost the same at this Mach number. An 
equation can be derived for the enstrophy from the Navier-Stokes equations." 1 It takes the following form by 
neglecting the dilatation and the viscosity fluctuations. 
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Here /^Production of enstrophy, r ffl =Diffusion of enstrophy by the turbulence, and £> ra =Diffusion of enstrophy by 
viscosity. S co is an important term that is responsible for the generation of vorticity fluctuations by the vortex 
stretching. It consists of three components: (1) the first term, S Mj i, consists of triple correlations among vorticity 
fluctuations and strain rate fluctuations, (2) the second term, S (!)> 2 , consists of correlations between the vorticity 
fluctuations and the stretching by the mean strain rate, and (3) the third term, S M> 3 , consists of correlations between 
the vorticity and the strain rate fluctuations. 


Figure 6 shows the root mean-square values of the vorticity fluctuations, co x>nm , co y rms , and co Zjrms for the channel 
flow. The maximum vorticity fluctuations for the streamwise (Ox.rms! an d spanwise (D z ,rms components occur at the 
wall. At the wall, for an incompressible flow 
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Hence the vorticity stretching terms, S m , the production term, the turbulent diffusion term T„„ and the advection 
terms are all zero at the wall. The dissipation equation becomes a balance between the molecular diffusion D m and 
the destruction term s m . Figure 7(a) and (b) show the different terms in the enstrophy equation. Instead of plotting 
the enstrophy, we plotted the equivalent kinetic energy dissipation quantity, 2 p*( enstrophy). The results show that 
near the wall, y + < 2, the viscous diffusion and viscous dissipation are large and balance each other. Within the 
region 2 < y + < 30, all the vortex stretching terms, Si, S 3 , and S 3 have the same magnitude contrary to the analysis 
by Tennekes and Lumley. In the core region y + > 30, the vortex stretching term Si becomes dominant and balances 
the viscous dissipation as inferred by the theory. These results agree with those presented in Mansour et al. 30 In Fig. 
7(b) we show the production term, the sum of the vortex stretching terms (S1+S2+S3) and the sum of the viscous 
dissipation and diffusion terms (s+T+D). The variations are much smoother than when considering them 
separately. It may be better to model the sums instead of modeling the individual terms. 

The standard model transport equation for the dissipation takes the following form 5 
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In RANS model development, 32 the terms in the dissipation equations are typically modeled as 
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Using the DNS data, Figs. 8(a), (b) and (c) compare these three modeled terms with “exact” terms computed 
directly from the enstrophy equation. Figure 8(a) shows that the shapes look similar, but the maximum amplitude 
obtained from the DNS data is about 50% higher than that predicted by the model. The agreement between the 
dissipation model and the DNS data is good except near the wall. This poor near-wall agreement is due to the fact 
that the model term becomes infinite near the wall. In many RANS models 3 ’ this deficiency is typically fixed by 
multiplying the model by a filtering function. The third term is very small, and the comparison is not good. This 
may be due to the inaccuracies when computing the term from the DNS data. 

B. Plane channel with constriction 


(a) Mean quantities 

The statistical quantities are obtained by averaging the solution in the spanwise direction and in time. The 
averaging in time was performed for about 15 flows through periods. This may be insufficient for separated flow 
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simulations. We are continuing the solution to run for more cycles. Figure 9 shows the mean streamwise velocity 
contours and the streamlines from the current DNS. We observe backflow and massive separation behind the first 
hill, as expected. Figures 10(a) and (b) show the skin friction C/ and the pressure coefficient c p along the lower and 
the upper surfaces. We also include the DNS results 17 obtained for the incompressible flow. In the current 
compressible computations at M 0 =0.2, the flow separates atx sep = 0.21 and reattaches at x reatt = 5.1. The separation 
point location is close to the DNS 17 and the LES simulations 18 predicted value of 0.21. The current predicted 
reattachment point location is about 5% upstream of that predicted by the DNS 17 and the LES 18 simulations, where 
the reattachment occurred at x reatt =5A and 5.3, respectively. We do not know the reason for these differences at this 
time. The agreement of the skin friction coefficient and the pressure coefficient from the present simulations with 
the DNS 17 are quite good. The skin friction coefficient along the upper surface is higher compared to the 
incompressible DNS 17 results. The lower surface skin friction coefficient becomes negative after the separation and 
reaches a minimum value near the middle of the separation bubble. Beyond this minimum point it slowly increases 
and takes a dip at the foot of the hill. The skin friction increases to larger values over the windward part of the hill 
and decreases at the crest. The pressure first increases and remains flat up to x ~ 2, and increases strongly up to the 
foot of the windward part of the hill, x ~ 7.5, before it decreases to the inflow value. 

Figures ll(a-e) show the average statistical quantities: (a) turbulent production, (b) turbulent dissipation, (c) 
ratio of turbulent production to dissipation, (d) turbulent kinetic energy, and (e) turbulent shear stress. Some 
oscillations are present due to insufficient sampling time. These quantities play important roles in RANS modeling. 
The production is maximum near the start of the separated shear layer and is confined to a region along the 
separation line. Beyond the reattachment point, the production is concentrated in the middle of the channel at a 
height ofy ~ 1. There is negative production near the reattachmcnt region and along the shoulder of the hill. The 
maximum dissipation occurs along the wall on the windward side of the hill. Outside the wall region, dissipation is 
confined to the shear layer and to its proximity. The ratio of production to dissipation is shown in Fig. 11(c). This 
ratio takes a peak value of 7 near the start of the separated shear layer. It decreases to a value around 2 in a region 
at a height ofy ~ 1. Recall that this ratio is about 1.8 in the buffer region of a turbulent channel flow. Figure 11(d) 
shows the turbulent kinetic energy distribution in the separated turbulent boundary layer. The kinetic energy is 
concentrated along the shear layer and it peaks near x ~ 3. After the reattachment point the energy is confined to 
region near y ~1. This agrees with the earlier observations from the turbulent production contours results. After the 
reattachment point, turbulence is confined to a region away from and parallel to the wall at a height ofy ~ 1. The 
turbulent shear stress contours in Fig. 11(e) also convey the same picture that turbulent shear stress is concentrated 
near the start of the shear layer and near they ~ 1 region at the foot of the hill. Figures 12(a-e) show the comparison 
of the mean streamwise velocity profiles with those obtained from the incompressible DNS 17 . The agreement is 
quite good except nearx = 4, where we obtained a smaller recirculation zone compared to Ref. 17. 

(b) Enstrophy computations 

Figure 13 shows the variation of different terms in the enstrophy equation at streamwise stations x = 0, 1, 2, 4, 
and 8. Instead of plotting the enstrophy, we plotted the equivalent kinetic energy dissipation quantity, 
2p*( enstrophy). There are some oscillations appearing in the profiles, due to lack of sampling and lack of accuracy 
in estimating the higher order correlations that appear in the enstrophy equations. We will concentrate on the 
profiles near the lower wall and in the separated shear layer region. At the wall, similar to the plane channel flow 
case, the dissipation and diffusion are much larger than other terms and balance each other. Although the figures do 
not show it because of the scales used, the magnitudes of these terms are on the order of 2 to 3 at all the stations. 
However, these terms decrease steeply to values on the order of 0.005 within a short height from the wall as seen in 
the figures. When we move away from the wall into the separated region, x = 1, 2 and 4, the viscous dissipation 
term, D t0 , and the vortex stretching term, Si u , gradually increase and peak locally in the shear layer. These terms 
peak close toy = 1.02, 0.95 and 0.68 at x = 1, 2, and 4, respectively. The vortex stretching component, S 2 M , exists at 
x=l, but its amplitude decreases downstream. Hence in the separated shear layer region the vortex stretching term, 
Sic, and the viscous dissipation term, D,,,, are the dominant terms in the enstrophy equation. This concurs with 
Tennekes and Lumley’s analysis that vortex stretching transfers the energy to the small scales and viscosity 
dissipates this energy. It is also interesting to observe that above the shear layer region the profiles look similar to 
those observed for the plane channel flow. These observations suggest that in free-shear flows such as a mixing 
layer, wakes, and jets, the triple correlation term will approximately balance the viscous dissipation term 
throughout the whole domain. 

(c) Computations using two-equations and Reynolds stress modeling 

To demonstrate the failings of RANS models for this type of separated flow, we also performed computations 
for this case using the two-equation k-co SST 33 model and the full Reynolds Stress Model (RSM) of Wilcox 5 . The 
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computations are carried out using the NASA code CFL3D. 34 Periodic boundary conditions are used in the 
streamwise direction. Figures 14(a) and (b) show the Mach number contours and the streamline patterns obtained 
with the SST and RSM models, respectively. Figure 15 shows the variation of the skin friction coefficient on the 
lower surface obtained from the RANS and RSM calculations. We also include the results obtained from the 
present DNS calculations and from Ref. 17. Although the separation location is predicted reasonably well by the 
models, the reattachment location is predicted too far aft, and the shape and the magnitudes of the skin friction 
coefficient compared to the DNS calculations are poor. The predicted reattachment points are located at x reatt ~ 7.8 
and 7.2 for the SST and RSM models, respectively. It is well-established that most RANS models typically under- 
predict the turbulent shear stress magnitudes in smooth-body separated shear layers, resulting in too little mixing 
and consequently delayed flow reattachment and boundary layer recovery. This is the case here as well. Rumsey 35 
and Rumsey and Jeyapaul 36 explored ideas that sensitize RANS models to regions of high P/s, resulting in earlier 
separated flow reattachment, but the fixes to date have been ad hoc. Further investigations along these lines are 
needed. 


V. Conclusions 

We investigated the turbulent flows through a plane channel flow and a channel flow with a constriction (2-D 
hill) using DNS and RANS calculations. We solved the N-S equations in the DNS using a higher order kinetic 
energy preserving central scheme and a fifth order accurate upwind biased WENO scheme for the space 
discretization. The kinetic energy preserving scheme did not work for the channel flow with the constriction and 
the solutions were obtained with the WENO scheme. We evaluated the magnitude of different terms that appear in 
the enstrophy equation. 

The mean flow profiles and the statistical quantities for the channel flow agreed with the available DNS results 
of Kim et al. 3 The analysis of the terms in the enstrophy equation for the channel flow revealed that the averaged 
vorticity dynamics could be divided into three regions. One is at the wall where the viscous diffusion and the 
viscous dissipation become large and balance each other. All the other terms are nearly zero at the wall. Away from 
the wall and below the buffer region, all three vortex stretching terms have a similar magnitude and their sum is 
balanced by the dissipation term. Beyond the buffer region, the triple correlation term among the vorticity 
fluctuations and the velocity strain perturbations dominates the kinematic part of the equation. This term is 
balanced by the viscous dissipation term. This concurs with the Tennekes and Lumley’s hypotheses that the energy 
from the large scales are transferred to the small scales by the vortex stretching mechanism and this energy is 
dissipated at the small scales by the viscosity. 

The simulation for the separated flow over a 2D-hill predicted a reverse flow and a large separation behind the 
hill. The flow separates slightly downstream of the hillcrest and reattaches in between the hills. The separation 
point location agrees with other DNS and LES simulation results. The computed reattachment point is about 5% 
upstream of those predicted by other simulations. We are pursuing the simulation with a finer grid to detect the 
reasons for this discrepancy. The turbulent statistical quantities revealed that the turbulent production reaches 
maximum values near the start of the separated shear layer. The turbulent dissipation always reaches maximum 
values at the wall. The ratio of production/dissipation reaches values on the order of 7 in the shear layer near the 
separation region. This is much larger compared to the plane channel flow, where it is about 1.8. The turbulent 
kinetic energy and Reynolds stress contours showed that high turbulence exists in a region parallel to the wall at 
approximately the hill height. This suggests that the shear layer formed at the separation point persists parallel to 
the wall up to the second hill. 

The enstrophy equation evaluation revealed that the vorticity dynamics is different in separated flows compared 
to that in equilibrium flows. The viscous dissipation and diffusion terms at the walls are about 100 times larger than 
in the outer part of the flow. These terms steeply decreases to very small values within a short height above the 
wall. Away from the wall in the separated region, the vortex stretching term and the dissipation term gradually 
increase, peak locally near the shear layer, then decrease outside of the shear layer. Beyond the shear layer, the 
profiles are similar to those observed in the plane channel flow. 

The computations with the RANS equations agreed with the previous findings that typical two-equation and 
Reynolds stress models predict a delayed reattachment point of the separation bubble. The question of how to 
extend or apply the findings from the analysis of the enstrophy and turbulent dissipation equations to turbulence 
modelling will be pursued in our future work. 
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(a) Plane channel 



(b) Channel with a 2-D hill 



Figure 1. Geometries and the coordinate systems used. 




Figure 2. Channel flow mean streamwise velocity (a) non-dimensionalized by the bulk velocity and the 
channel half-width, (b) non-dimensionalized by viscous units. 




Figure 3. Channel flow root-mean-square of the velocity fluctuations and the shear stress normalized by the 

wall shear velocity. 
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Figure 4. Channel flow balance of different terms in the turbulent kinetic energy equation. 


(a) uu (b) vv 




(c) ww 


(d) uv 




Figure 5. Channel flow balance of different terms in the Reynolds stress equations for (a) u’u’, (b) v’v’, (c) 

w’w’, and (d) u’v\ 
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Figure 6. Channel flow root-mean-square values for the vorticity fluctuations. 


(a) 


(b) 




Figure 7. Channel flow balance of different terms in the enstrophy equation, 
(a) (b) (c) 



Figure 8. Comparison of terms in the modeled dissipation with the corresponding terms in the enstrophy 
equation. Symbols - modeled terms, solid lines - from the enstrophy equation; (a) S £2 - S e3 , (b) D £ - S £ i; (c) T £ . 
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Figure 9. Hill flow contours of the mean streamwise velocity and the streamlines, Re=2800. 


(a) Skin friction (b) Pressure 




Figure 10. Hill flow variation of mean (a) skin friction coefficient, and (b) pressure coefficient along the 

lower and upper walls. 
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Figure 11. Hill flow contours of the turbulent statistical quantities (a) Production, (b) Dissipation, (c) 
Production/Dissipation, (d) Kinetic energy, and (e) Shear stress. 
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(a) x=0 (b) x— 1 (c) x=2 (d) x=4 (e) x=8 



Figure 12. Hill flow mean streamwise velocity profiles at different stations. 

(a) x = 0 (b) x = 1 (c) x = 2 





(d) x = 4 (e) x = 8 




Figure 13. Hill flow balance of different terms in the enstrophy equation. 
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Figure 14. Hill flow streamwise velocity contours and the streamline patterns obtained from RANS 

calculations, (a) SST, (b) RSM. 



Figure 15. Hill flow variation of the skin friction along the lower surface obtained from (1) present DNS, (2) 
from the incompressible DNS 17 , (3) SST, and (4) RSM models. 
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